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We extend previous studies of the BCS canonical approach for the attractive Hubbard model. A 
derivation of the BCS formulation is presented for both the Hubbard and a simpler reduced Hamil- 
tonian. Using direct diagonalization, exact one and two dimensional solutions for both Hamiltonians 
are compared to BCS variational calculations. Approximate and exact ground state energies and 
energy gaps are computed for different electron number systems as well as correlation observables 
not previously predicted. Reproducing published one dimensional findings, the BCS is an excellent 
approximation for the Hubbard ground state energy but not energy gap, a finding that remains 
true in two dimensions. Propagators and correlators are found more sensitive to wavefunctions and 
appreciable differences are computed with the Hubbard model exhibiting a weaker degree of super- 
conductivity than the BCS. However for the reduced Hamiltonian model the BCS is an excellent 
approximation for all observables in both one and two dimensions. 

PACS numbers: PACS number(s): 74.20.Fg, 71.24. +q, 71.10Fd, 71.10Li 



I. INTRODUCTION 

Since the advent of the theory of superconductivity, de- 
veloped by Bardeen, Cooper and Schrieffer [BCS}i and 
their contemporaries, Bogoliubov- and ValatinSi, count- 
less investigations with these formalisms have been per- 
formed in a variety of fields. This widespread application 
of the BCS approach is due not only to its power and 
utility but also because superior approximate or exact 
solutions are generally not available. This is especially 
true in nucleas^ and hadronic structure^ where physical 
insight has been limited by uncertainties due to approxi- 
mations. Hence, further assessing the validity of the BCS 
approximation is an important endeavor with broad at- 
tending interest. In particular, there have been several 
studies£i2i£ of the BCS approximation within the Hub- 
bard model for which exact solutions, such as those pro- 
vided by the Bethe Ansataifi, are known but only in one 
dimension. The purpose of this paper is to extend these 
investigations in several new directions: 1) assess the 
BCS approximation for not only the Hubbard but also 
a reduced Hubbard-like Hamiltonian; 2) in addition to 
energies and energy gaps also examine propagators and 
correlation observables; 3) evaluate the BCS accuracy in 
two dimensions. 

In discussing the Hubbard model it is natural to re- 
strict ourselves to a fixed filling because the Hamilto- 
nian conserves electron number. Therefore we utilize 
the canonical BCS formulation with a fixed number of 
electrons. This is superior to the grand canonical ap- 
proach which only constrains the average electron num- 
ber and also produces a higher ground state energy. A 
previous pioneering study^ notably detailed the excel- 
lent ground state energy agreement between the exact 
Hubbard model and the approximate canonical BCS so- 
lution. However that work also determined, depending 
upon electron density, that the Hubbard and BCS en- 
ergy gaps do deviate, again in one dimension. Related, an 
odd filling BCS application^ in one dimension also found 



somewhat less precise results. Because the BCS approx- 
imation is conjectured to improve in higher dimension 6 , 
it is important to perform a numerical assessment in two 
dimensions. 



When computing with the BCS ground state many of 
the terms in the Hubbard Hamiltonian, H, become ir- 
relevant since only components which are either diago- 
nal in momentum space or create and destroy momen- 
tum modes in pairs contribute. Truncation to these rem- 
nant terms defines the reduced Hamiltonian, H r , which 
is specified in the next section where the BCS Hubbard 
ground state energy, EbcSj is found degenerate with 
the BCS reduced Hamiltonian value. Because H con- 
tains further terms describing additional dynamics we a 
priori expect Ebcs to better approximate the exact re- 
duced Hamiltonian ground state energy, E r , than the 
exact Hubbard ground state energy, E. Indeed, as we 
detail below Ebcs is an order of magnitude closer to 
E r than E even though the variational treatment is very 
good for both. Consequently a more complete BCS as- 
sessment, especially for systems modeled by H r , should 
address other observables which, in contrast to energies, 
are more sensitive to wavefunction details. Accordingly, 
we also calculate electron pair propagators, wavefunction 
inner products and s, p and d-wave correlators which are 
important for establishing superconductivity. 



This paper is divided into four sections. In section II 
the Hubbard and reduced Hamiltonians are introduced 
and a new formulation of the canonical BCS method is 
derived. This section also details formulas for energy 
gaps, densities, propagators and s, p and d-wave cor- 
relators. Section III compares and discusses numerical 
results for a variety of electron level fillings and different 
number of lattice sites. Finally, conclusions are summa- 
rized in section IV. 
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II. THEORETICAL FORMULATION 
A. Model Hamiltonians 

We follow the literature notation, especially Refs.^and 
and recall the full Hubbard Hamiltonian in momentum 



space 



H 



u sr t t 

pp'q 



(1) 



The reduced Hamiltonian retains the non-diagonal (p ^ 
p') terms with q = and diagonal (p = p') terms 

H r = £p a U a P<T + 77 E a P T a -pi a -p'l a P'T 



pa 



V 

u 



pp'p^p' 



E«: 



X a -p+qi a pT ■ ( 2 ) 



pq 

In the above equations a pa .(a pa ) creates (annihilates) an 
electron with momentum p and spin a = f or |. In one 
dimension the electron kinetic energy is e p = — 2icos(p), 
while in two dimensions, e p = — 2t(cos(p x ) + cos(p y )), 
with hopping rate t (set to 1 in this work) between near- 
est neighbors. The momenta are restricted to multiples 
of ^r, where L is the linear extent of the system. The 
number of sites V is L in one dimension or L 2 in two 
dimensions. Note that the last term in Eq. J2J can 
be expressed in terms of a sum over the dummy index 
1 = p + q yielding ^N e ^N e i, the product of spin up, 
iVef, and down, N e ±, electron number operators for a 
system having V total lattice sites. Since we restrict this 
discussion to the attractive Hubbard model the potential 
coupling parameter U is negative. The reduced Hamilto- 
nian, H r , has been chosen to only contain elements of H 
which have non-zero expectation values in the BCS state 
defined in the next section. 



B. BCS Approximation 



I* 2 -) = e n^°u«- 

s(l/:) i=l 



10} • 



(5) 



Here and below the cryptic summation notation s{v :) is 
used to represent the sum over all unordered sets of mo- 
mentum coordinates, k 1; k 2 , . . . k^, with the restriction 
that all of the momenta are different. 

For a physical system with N e electrons the BCS 
ground state energy Ebcs, is calculated by minimizing 
the expectation value E 2v 



(6) 



i.e. Ebcs — mm E 2v where the minimization is achieved 
by varying the g^ as detailed below. 

Before proceeding several points are in order. First, in 
the grand canonical BCS formulation the ground state 
energy is obtained by minimizing ^bcs^ekjs)'* ■ Using 
Eq. Js3 one can show that the grand canonical is greater 
than the minimum canonical energy. Second, combining 
Eqs. 11121')!) directly yields the degenerate result 



2ul 



(*2»\H r \* 2v ) 
(^2u\^2u) 



(7) 



verifying, as stated above, that the BCS ground state en- 
ergy is the same for both Hamiltonians. This is also true 
for the grand canonical BCS energy. Third, the succes- 
sive operation of H r on |0) only generates states having 
electrons paired just like the state |BCS) with opposite 
momentum and spin, called Cooper pairs. In general, 
the exact Hubbard ground state will contain additional 
electron correlations between such pairs and the BCS ap- 
proach assumes these correlations can be ignored. 
Now, we compute the norm of the fock state 



s{w.) i=l 

and the unnormalized kinetic energy 



(8) 



We adopt the formulation of Ref. |6| but present an al- 
ternate BCS derivation. For a given Hamiltonian, the 
BCS approximation entails using the unnormalized vari- 
ational wavefunction 



\BCS) =H(1+: 9k4 T «-ki)|0) , 



(3) 



to minimize the ground state energy. Here are the vari- 
ational parameters and |0) represents the vacuum state 
with no electrons. Note that the above Hamiltonians 
conserve particle number, but the state \BCS) includes 
a mixture of fock components with different particle num- 
ber. Because we are interested in states with fixed fill- 
ings we organize these components by electron number 
iV e = 2v with v = 0,1,2,3... 



\BCS) =Y,\^2 V ) 



(4) 



i^2u\J2 e P a U a p^^) = 2 H £ p5p e n 

P<? P s(v— l:p) i 



The modified summation index s(n : pi, P2 ■ ■ •) now de- 
notes the sum over all unordered sets of momentum coor- 
dinates, ki, k 2 , . . . k„, with the stronger restriction that 
all momenta in the extended set, {ki, . . . k„, Pi, p 2 , ■ ■ ■}, 
are different. For the potential energy we first calculate 



Y a - P l a Pl 1*2*) 



p s(w.) i=l 

= E^p E n^4 iT a -kai°) • ( 9 ) 

P s(v— l:p) i=l 
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Then the unnormalized potential energy matrix element 
is 

^l|(E a PT a -Pl E a-p'iVT+tf TAT 1 W„) 



17 

r 



p'^p 



p,p' s(i/-l:p,p') i=l 



Note well that the sum over the set s(v— 1 : p, p') ensures 
p ^ p'. Now define 

n 

A-( P1 ,... Pro )= 2 n^, (io) 

s(n:pi,p 2 ,...p m ) »=1 



+^E^p'^-i( p/ ' p ) 

2C7 



p' 



~y E fp'5p"5 P ^-2(P'>P">P) 



P':P" 



^gpAl^p^ - — ) . 



Solving for g p yields 



gp = T7^-^29p'^1-i(p',P) , 



VV, 



p' 



with the denominator, 



(16) 



(17) 



and 



B™(pi,...p m ) 



= A™( Pl ,... Pri 

A?, 



(11) 



which both vanish if any two momenta are equal, to ex- 
press E 2v as 



Ve = Al_ x { V ){E^- — -2e p ) 



V 



' 2 E £ p'.9p'^-2(p':P) 



E 9P'9 P "^- 2 (P'.P".P) • (18) 



E 2v -v 2 - 

2.U y 



(12) 



2E P e P .9p^-i(p) + |E P ,p- g P g P ' 4Li(P,P') 
2 E £ p5p^-i(p) + ^E» B -i(P'P')- 



Next, we evaluate the BCS variational equation for a 
state with N„ = 2v electrons 



dE 2v 



o , 



(13) 



by first calculating 

6>A™( Pl ,...p m ) 



e n £ 

s(n-l:p 1: p 2 ,...p m ,p) i=l 

2 gp A™+ 1 1 ( Pl ,..., Pm ,p) . (14) 



Then we use the product rule for differentiating, 



dB™(p u ...p m ) 

dg P 

dA^(p 1 ,...p m ) 



(15) 



/CiP:.-..p, 



dAl 



dg P ^'dg p 

2 5p «_+ 1 (p 1 , . . . p m ,p) - B™(Pl, ■ • • PmMi-l(p)), 



to obtain 
= A 



9i?2iy 



<9g F 



= 4epff p Ai_ 1 ( P ) +4^e p ^^5p^_ 2 (p / , P ) 
p' 



Equation (|17J1 is equivalent to Eq. (34) of Ref. (note 
their equation has a misprint involving a missing g p ). 

In the case of odd filling, N e — 2^ + 1, the variational 
state has a single specific unpaired, say up, electron with 
momentum q, chosen to minimize the energy: 



|*a«+i) = E ( II^<T at -ka KtI°> ' ( 19 ) 

s(>:q) M=l ' 

Operating with the Hamiltonian on this state shifts the 
energy by e q | + yN± and repeating the above analysis 
yields 

.U 

E 2v +i - e q - ^(f + 1) — 

2 E P e P gp^-i(p ; q) + fE P ,p- gpgp'^-iCp.p'.q) 



^i(q) 



u 



2E £ p5p C ^i(p)+ 7/ E^P'^-i(P.P') . (20) 



p>p' 



where 



W (Pi, ■ ■ -Pm) = 7773 > ( 2i ) 

is equivalent to i?™(pi, . . . p m ) defined on a space with 
the q mode removed. Finally, the odd electron system 
variational equation reduces to 

9p = y^- o E 9p ,A l-i (p'. p. <0 > ( 22 ) 

p' 

where 
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V = ^(Rq)^ - V y ' - 2e p - e q ) 

- 2 ^ £ p'5p'^- 2 (p' I p I q) 
p' 

"7 E .9p'.9 P "^- 2 (p',P",P,q) • (23) 
p',p" 

C. Exact diagonalization 

Using PARPACK, a parallel implementation of the 
Lanczos algorithm, we exactly calculated the Hubbard 
and reduced Hamiltonian ground states as a function of 
the electron density. The dimension of the Hilbert space 
for a 16 site Hubbard model is 4 16 w 4 x 10 9 . Since the 
Hamiltonians conserve momentum and particle number 
we may restrict to subspaces with fixed momentum and 
filling. However, the dimension for half- filling at mo- 
mentum is still greater than 10 7 . There are about 2000 
non-zero terms for each state, so the Hubbard Hamilto- 
nian matrix cannot be stored in memory even in a sparse 
format. Instead, matrix elements were generated at each 
step in the Lanczos algorithm. The diagonalization took 
about 3 hours using 30 1.2 Ghz Athlon processors in 
parallel. 

For even N e = 2is, we restrict to the zero momentum 
subspace which yields the minimum energy eigenvector. 
For odd N e , the momentum with lowest ground state 
energy is dependent on the coupling U. For U = it is 
obvious that the lowest energy will have the momentum 
of the v + 1th lowest energy mode. For high U, with less 
than half filling, the lowest energy will have momentum, 
but for the value U = — 10 used throughout this paper, 
the momentum is the same as for U = 0. 

The total momentum of the ground state for odd N e 
corresponds to occupation of a specific unpaired momen- 
tum mode. For the Hubbard model, occupation of this 
mode is not fixed by the action of the Hamiltonian but, 
for the reduced model, all relevant basis states have the 
same unpaired mode. In this regard, the reduced model 
ground state is the same as the BCS state which has the 
single unpaired mode q (Eq. IT§|| . 

D. Observables and theoretical constructs 

To properly assess any approximation scheme, as well 
as investigate the robustness of a particular dynamic 
model, it is important to calculate and compare a va- 
riety of observables and theoretical constructs. 

In addition to the ground state energy, E, and energy 
density, E/V, another important observable is the single 
particle energy gap, Ajv e , 



= ~[E Ne+ i - 2E Ne + E Ne -!} . (24) 

In the grand canonical formulation, a quasi-particle exci- 
tation is produced by the action of , or equivalently, 
a_ p j on the state |BCS)— with excitation energy 

_ (BCS|Qp T gq^ T |BCS) (BCS|ff|BCS) 
p ~ (BCS\BCS) (BCS\BCS) 

= V A bcs+4 ■ (25) 

Here Abcs is the BCS gap and e p is dependent on the 
momentum p. The minimum excitation is non-zero since 
E p > Abcs i thus accounting for the name "gap" . In the 
canonical BCS formulation Op| x f , 2i/ is no longer propor- 
tional to a P i^>2v Therefore, An s averages the effects of 
adding and subtracting one electron. 

Because wavefunctions are more sensitive than ener- 
gies to approximation schemes, we also address other 
operator wavefunction overlaps such as propagators 
and correlators 9 . Before introducing theses constructs 
we first specify the generalized one, p aT (r,s), and 
two, p . ./ rT ;(r, r', s, s'), particle density matrices for any 
model state with unit norm 

p ffT (r,s) = {y\clc ST \M>} = ( 4 CT c ST ) (26) 

/Wi-r'(r,r',s,s') = ( 4 ct 4' CT ' c stCs't' } . (27) 

Here the indices r, s represent possible electron lattice 
sites. The position operators, c rcr and cJ CT , are related to 
the momentum operators by 

<V =^E P ^ pr a P (28) 

Between states with the same filling pjj(r,s) = 
Pxi(r, s) = 0, while p^(r,s) and p^(r,s) are the one 
electron propagators for spin up or down, respectively. 

Next we construct a variety of correlators as the dif- 
ference between a specific two particle density and the 
product of one particle densities. For example, the on 
site s-wave correlator, C s (r,s), is 

C s (r,s) = p TUT (r,r,s,s) - p TT (r, s)p u (r, s) , (30) 

which is the correlated electron pair propagator minus 
the product of two (independent) one electron propaga- 
tors. The extended s-wave correlator, C es (r,s), is 

C es (r,s) = ^2w es (t)w es (t')(p n ^(r,r + t,s,s + t') 
t,t' ^ 

-PTTMVu( r + M + t')) , (31) 
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which involves a sum, t = (t x ,t y ), t' = {t' x ,t' ), over 
the four planar sites adjacent to lattice site r (s) with 
weights w es (t), w es (t') given by w es (t) = 1 for t = 
(1,0), (0,1), (-1,0) or (0,-1) and iu es (t) = otherwise. 

We also list in Table I several other correlators, C sym , 
with different symmetry weights w sym (t) which are im- 
portant for investigating superconductivity. See Ref. 9 
for illustrations and further discussion. 



TABLE I: Correlators 



Correlator 


C sym 


non-zero w S ym(t) terms 


s-wave 


C 8 


10,(0,0) = 1 


extended 
s-wave 


Ces 


tu es (l,0) =w es (-l,0) = 1 

10 M (0,1) =Wea(0,-l) = 1 


x direction 
p-wave 




Ufa (1,0) = 1 
w P J-l,0) = -l 


d-wave 


Cd ;c 2_ y 2 


W d„2_ y 2 (1, 0) = W d^2_ y 2 ("I, 0) = 1 
W <* x 2_„2 (0, 1) = Wd x2 _ H2 (0, -1) = -1 


diagonal 
<i-wave 




^(1,1) =w dxy (-1,-1) = 1 
«Jd„(l,-l)=to*..(-l,l) = -l 



Since our system has translational symmetry, the cor- 
relator, C sym (r, r + d), and propagator, p aa (r, r + d), do 
not depend on r. Therefore we can write 

C sym (d) = C sym (r, r + d) = — ^ C sym (s, s + d) (32) 

S 

Paa(d) = Paa(r, T + d) = — ^ P<?<? ( S > S + d ) • ( 33 ) 

s 

As discussed in Ref. H an effective measure of super- 
conductivity is the concept of Off-Diagonal Long Range 
Order (ODLRO) which can be quantitatively character- 
ized by the correlator or ODLRO function, C sym (d). 
The criterion for demonstrating superconductivity is that 

C sym 

(d) remains finite, positive and does not decrease to 
zero at large r. 

Using Eqs. (|28l 129(1 . the propagator can be written in 
momentum space as 

PTt( d ) = y^2( <4| c r+dT ) 
r 

= T^E^P'^Kt^T > , ( 34 ) 

where 

P( P ,q,d) = ^y(PT-q-(r+d)) 

r 

= U<5 pq e- iq d . (35) 

Thus 

PTT(d) = ^E e ~ !P ' d <°PT a PT ) - ( 36 ) 
p 



and evaluating C syrn (d) in momentum space yields 

(t)w sym (t r ) (37) 

r,t,t' 

p mT (r,r + t,r + d + t',r + d) 

= Va E F (P' P'' P"> P'"' d )< a pT a P'l a p"i a P'"T ) (38) 
p,p',p",p"' 

P(P, p', p",p"', d) - ^ p , p , +p ^_ p , e - 4d (p"+P'") (39) 

t t' 
In particular for d-wave symmetry 

P(p,p',p",p'",d) = W8 p , p n +p ,„_ p ,e- lA ^"+v"') 
(cos(^) - cos(p;)) (cos(p£) - cos(^)) . (40) 

The other correlators were calculated similarly. 



III. NUMERICAL RESULTS AND DISCUSSION 

A. One dimensional systems 

We first compute for a one dimensional lattice the ap- 
proximate BCS and exact Hubbard and reduced Hamil- 
tonians ground state energies as a function of level filling. 
In Fig. ^the Hubbard (boxes), E/V, reduced (circles), 
E r /V, and BCS (crosses), Ebcs/V, ground state energy 
densities are compared as a function of electron filling, 
N e , for t = 1, V = 16 sites and U = -10. The BCS and 
exact Hubbard data reproduce the data in Refs.0 and0- 
As they reported, the agreement is very good, especially 
for the case of even filling. The exact reduced data is 
indistinguishable from the BCS data in the figure. Over 
99% of the difference between the BCS and the exact 
Hubbard is due to the difference between H r and H and 
only 1% is due to quantum correlations between different 
BCS pairs which are excluded by the variational wave- 
function. Similar results have been obtained for other 
potential strengths and different number of lattice sites 
(not shown). 

The Hubbard (boxes), reduced (circles) and BCS 
(crosses) energy gaps, A Ne , are compared in Fig. [21 As 
with the ground state energy, the BCS gap is the same 
for either the Hubbard or reduced Hamiltonian. In con- 
trast to energy density comparisons, the BCS-Hubbard 
difference is more pronounced especially for half-filling, 
N e = V. The BCS is still a good treatment for the re- 
duced Hamiltonian. 



B. Two dimensional systems 

Having reproduced and extended one dimensional 
studies, we now address two dimensional systems. Figure 
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FIG. 1: One dimensional ground state energy densities as a 
function of electron filling for the Hubbard (boxes) and re- 
duced (circles) Hamiltonians. Note the good agreement with 
the BCS results (crosses) . 



FIG. 3: Two dimensional ground state energy densities as a 
function of electron filling for the Hubbard (boxes) and re- 
duced (circles) Hamiltonians. Note the good agreement with 
the BCS results (crosses). 
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FIG. 2: One dimensional energy gaps as a function of electron 
filling for the Hubbard (boxes) and reduced (circles) Hamil- 
tonians. The BCS (crosses) is least accurate near half-filling. 

3 depicts the Hubbard, reduced and BCS energy densi- 
ties for a two dimensional lattice having y = 4x4 = 16 
sites. 

Again the BCS energies are close to the exact Hubbard 
eigenvalues and indistinguishable from the reduced eigen- 
values. The corresponding energy gap comparisons are 
displayed in Fig. 4. Contrary to conjecture, the differ- 
ence between the exact and BCS energy gaps remains 
similar to the one dimensional case. 

We further note that the Hubbard and reduced energy 
eigenvalues are generally in reasonably good agreement 
for all electron densities in both one and two dimensions. 
Accordingly, it would appear that the reduced Hamilto- 
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FIG. 4: Two dimensional energy gaps as a function of elec- 
tron filling for the exact Hubbard (boxes) and reduced (cir- 
cles) Hamiltonians. Note the disagreement between the BCS 
(crosses) and the Hubbard energy gaps, similar to the one 
dimensional result (compare to Fig. 2) . 

nian incorporates the dominant physical dynamics. This 
is fortuitous because H r can be solved more easily than 
H, whether by the BCS approximation or by diagonaliza- 
tion, which in turn facilitates investigating more complex 
systems in multi-dimensions. 

C. Wavefunction overlap 

In general, a variational calculation reproduces the 
ground state energy much better than other observables. 
It is therefore of interest to examine the overlap of the 
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BCS and exact model wavef unctions. Figure [5] plots the 
inner product of the exact, one dimensional Hubbard 
eigenvector with the BCS state as a function of elec- 
tron filling. Notice the significant deviation from unity 
especially with increasing density. As before, the BCS 
and Hubbard results have the largest difference for sys- 
tems near half-filling and odd number of electrons. In 
two dimensions the overlap behavior is similar but now 
the values are significantly higher (see Fig. |HJ, consis- 
tent with the conjecture that the BCS approximation 
improves with higher dimensionality. The inner product 
of the reduced model eigenvector with the \BCS) is .999 
in both 1 and 2 dimensions (not shown). 




10 12 14 16 



FIG. 5: Inner product of the Hubbard and BCS vectors as a 
function of filling N e (one dimension). 
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FIG. 6: Inner product of the Hubbard and BCS vectors as a 
function of filling N e (two dimensions). 



D. Densities, propagators and correlators 

Returning to physical observables, this section com- 
pares propagators and correlators for exact and approx- 
imate model wavefunctions. Unless stated otherwise, all 
calculations were performed for U = —10, lattice sizes 
V = 16 x 1 and V = 4 x 4, and N e = 15 and iV e = 16. 



Results for other fillings were also calculated and found 
to be qualitatively similar. 

Figures and |H1 compare the Hubbard, reduced and 
BCS electron propagators, p^(d), for different lattice 
ranges d and N e = 16 in one and two dimensions, respec- 
tively. The BCS and reduced propagators are essentially 
identical and also similar to the Hubbard result. At lower 
densities there are somewhat larger differences between 
the Hubbard and BCS (or reduced) propagators. 
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FIG. 7: One electron propagators as a function of distance 
(one dimension, Ne = 16). 
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FIG. 8: One electron propagators as a function of distance 
for (two dimensions, N e = 16). 

Figures El and ED compare the different model sin- 
gle electron propagators again in one and two dimen- 
sions, respectively, but now N e = 15. Notice that, just 
like the model energy comparisons, the Hubbard and re- 
duced propagators are more different for odd N e , but 
again the BCS accurately reproduces the reduced results. 
Although lower N e Hubbard-reduced comparisons yield 
even larger devations (not shown), overall there is no sig- 
nificant difference in agreement between the two propa- 
gators in both one and two dimensions. 
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FIG. 9: One electron propagators as a function of distance 
(one dimension, JV e = 15). 
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FIG. 10: One electron propagators as a function of distance 
(two dimensions, N e — 15). 



Figures ITT1 and IT21 similarly compare the Hubbard, re- 
duced and BCS electron s-wave correlators, C s (d), for 
N e = 16 in one and two dimensions, respectively. Now 
there are significant departures between the Hubbard and 
reduced results. Notice that the reduced correlator has a 
larger magnitude than the Hubbard value. This is consis- 
tent with the energy gap comparison, where the reduced 
gap was larger than the Hubbard gap. Also, but not 
shown, lower N e yields smaller Hubbard-reduced corre- 
lator differences. The BCS again very accurately repro- 
duces the reduced correlators. 

The Hubbard and reduced s-wave correlators are closer 
in two dimensions than in one dimension but the lim- 
ited range of d sampled warrants caution in this conclu- 
sion. By the Mermin- Wagner theorem^, the Hubbard 
model will not display long range correlations at T ^ 
in one or two dimensions. However, according to the 
Bethe ansatai^, at least for large \U\, the ground state 
energy of the one-dimensional system becomes part of 



a continuous spectrum as the system length increases. 
The theorem does not apply at T = under these cir- 
cumstances and therefore long distance correlations are 
possible. The limited data shown in figures ITT1 and IT21 are 
inconclusive on the question of long distance order. 

Conversely, the BCS and reduced ODLRO functions 
are flat for increasing d. The Mermin- Wagner theorem 
does not apply to the reduced Hamiltonian, even at T ^ 
0, because it has long range interactions. For the ground 
state of the reduced Hamiltonian, the s-wave correlator 
reduces to 
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(41) 



The second term decreases linearly with V for large d, 
but the first term is independent of d and approaches a 
constant as V increases. For U,N e ,V and \d\ 3> 1, this 
constant is Ne<y2 ^ N ^ ■ A precise limit can also be given 
for the BCS correlator at finite V. Resummation of the 
correlator yields terms independent of d, 



p^p' 



b=l 



-l) b ]T ffp ^_ 6 (p), 
p 



with additional terms decreasing rapidly to with in- 
creasing d. 
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FIG. 11: The s-wave correlators as a function of distance (one 
dimension, N e = 16). 

The Hubbard, reduced and BCS electron s-wave corre- 
lators, C s (d), for N e = 15 are compared in Figs. 1131 and 
1 141 for one and two dimensions, respectively. The differ- 
ences are more dramatic for odd fillings as the Hubbard 
correlator in one dimension falls rapidly to 0. The two 
dimensional Hubbard correlator also falls rapidly within 
the limited range in d. In contrast, the BCS and reduced 
models, with a single specific unpaired electron mode, 
still exhibit long range correlations. 
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FIG. 12: The s-wave correlators as a function of distance (two 
dimensions, N e = 16) 



FIG. 14: The s-wave correlators as a function of distance (two 
dimensions, N e = 15). 
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FIG. 13: The s-wave correlators as a function of distance (one 
dimension, N e = 15). 



FIG. 15: The p-wave correlators as a function of distance (one 
dimension, N e = 16). 



The p-wave correlators, C Px (d), are plotted and com- 
pared in Figs. [13 CHI El and dl for different fillings 
and dimensions. Both interesting and somewhat surpris- 
ing is that the differences between the BCS and reduced 
correlators are now discernable, especially in 2 dimen- 
sions. Since these two model wavefunctions are identical 
to within 1 part in 10 3 , this indicates that the p-wave cor- 
relator is a very sensitive wavefunction probe as this cor- 
relator's matrix elements are dominated by a very small 
portion of the state vector. This result also reflects the 
small p-wave correlator magnitude. None of the models 
exhibit strong long range order or significant p x correla- 
tions. Finally, the differences between Hubbard and re- 
duced p-wave correlators are dramatic for both one and 
two dimensions with even or odd filling. 

Finally, the d-wave correlators, Cd aj2 _ 2 > were calcu- 
lated in two dimensions (Figs. 1191 and I2U|) . As with the 
p-wave correlators, the magnitudes are small for non-zero 
d and the differences between the Hubbard and reduced 



correlators are dramatic. The differences between the 
reduced and BCS correlators are again present. 
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FIG. 16: The p-wave correlators as a function of distance (two 
dimensions, iV e = 16). 



FIG. 19: The d-wave correlators as a function of distance (two 
dimensions, N e — 16). 
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FIG. 17: The p-wave correlators as a function of distance (one 
dimensions, iV e = 15). 
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FIG. 20: The d-wave correlators as a function of distance (two 
dimensions, N e — 15). 



FIG. 18: The p-wave correlators as a function of distance (two 
dimensions, iV e = 15). 
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IV. CONCLUSIONS 

The major thrust of this work has been to further 
numerically assess the accuracy of the BCS approxima- 
tion. By applying this approximation to multi-electron 
systems governed by both the Hubbard and a simpler, 
reduced Hamiltonian, each exactly diagonalizable, it has 
been documented that the BCS provides an excellent de- 
scription of ground state energies in both one and two di- 
mensions. However, the BCS single particle energy gaps 
are less accurate, in both one and two dimensions, es- 
pecially near half filling. This undermines previous as- 
sertions that the BCS approximation improves with in- 
creasing dimension. Further, propagators, overlaps and 
correlators have been found to be quite sensitive to both 
interactions and approximations as the BCS and reduced 
model predictions are appreciably different than the Hub- 
bard results. The slow decrease of the Hubbard and re- 
duced Hamiltonian ODLRO propagators and correlator 
functions with increasing lattice separation is evidence 
for superconductivity. However correlations, and thus 
the degree of superconductivity, are much stronger in the 
reduced model. Finally, the BCS method extremely ac- 
curately reproduces all reduced model observables. 

In summary, our results challenge the validity of the 



BCS method as a good, comprehensive approximation to 
the attractive Hubbard model. While the Hubbard and 
BCS correlators both exhibit long-range order, their vari- 
ations with distance are quite different. Because agree- 
ment between correlators and inner-product overlaps im- 
proves to some extent with dimension, perhaps in three 
dimensions the BCS method may yet be a very accurate 
alternative to exact Hubbard model wavefunctions but 
further study is necessary. Finally, for physical systems 
where the interactions embodied in the reduced Hamil- 
tonian dominate, the BCS approximation will be a very 
accurate approximation scheme permitting significantly 
more ambitious numerical investigations. 
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